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Abstract 

Cosmological evolution of axionic string network is analyzed in terms of field- 
theoretic simulations in a box of 512 3 grids, which are the largest ever, using a new 
and more efficient identification scheme of global strings. The scaling parameter 
is found to be £ = 0.87 ± 0.14 in agreement with previous results. The energy 
spectrum is calculated precisely using a pseudo power spectrum estimator which 
significantly reduces the error in the mean reciprocal comoving momentum. The 
resultant constraint on the axion decay constant leads to f a < 3 x 10 11 GeV. We also 
discuss implications for the early Universe. 



1 Introduction 



The axion has rich implications for astrophysics and cosmology. In origin, the axion is 
a pseudo Nambu-Goldstone (NG) boson of the Peccei-Quinn (PQ) symmetry, which is 
introduced to the Standard Model as a solution to the Strong CP problem in quantum 
chromodynamics (QCD) pQ. Below the QCD scale Aqcd — 200MeV, the axion acquires a 
mass from the instanton effect [2j[3], 

fa 



m a = 6 fieV 



10 12 GeV 



where f a is the axion decay constant. Couplings to other particles are also suppressed by f a , 
which is constrained from various terrestrial experiments and astrophysical observations. 
The most stringent lower bound f a > 10 10 GeV comes from the duration time of neutrino 
burst from SN 1987A HE]. 

On the other hand, f a is also bounded from above based on cosmological considerations. 
While couplings of axions to other particles are very week, they can be copiously produced 
in high-energy environments in the early Universe. These axions contribute to the energy 
density of dark matter. If inflation [6], which is essential in current understanding of our 
Universe, takes place in the early epoch and the PQ symmetry is broken during and after 
inflation, the axion field becomes homogeneous over the observable Universe apart from 
quantum fluctuations, and undergoes coherent oscillation after the QCD phase transition 
with an initial value 0(f a ). Such coherently oscillating axion behaves as cold dark matter 
(CDM), and should not dominate over the observed energy thereof. From the cosmic 
microwave background (CMB) and other cosmological observations, the abundance of 
the CDM is precisely constrained, which gives a current bound on the decay constant 
fa as fa 10 11 GeV [7]. If the PQ symmetry is not restored after the inflation, the 
isocurvature fluctuation of axion leads to significant CMB anisotropies. The amplitude of 
axion isocurvature perturbation is proportional to the Hubble parameter in the inflationary 
epoch -f/i n f, and recent observations give a constraint H ini < 10 7 GeV for f a ~ 10 10 ~ 12 GeV 
and misalignment angle 6 a ~ 0(1) [5HTT]. Realizing such an inflation model is generally 
difficult and requires fine tuning of parameters. This calls up another possibility, where 
the PQ symmetry is restored during or after the inflation. This happens rather naturally 
in the supersymmetric inflation models where the effective mass of 0(H m {) is induced for 
the PQ scalar through supergravity effect and (if it is positive) restores the PQ symmetry. 
Other models to restore the symmetry during inflation have been proposed in [T2|[T5]. 

Since the PQ symmetry is a global U(l) symmetry, linear topological objects called 
axionic strings can form when it breaks spontaneously. When the spontaneous breaking 
of U(1)pq occurs in the Universe, a cosmological network of axionic strings is formed. 
The axionic string is a global string and an extended object with physical width d str mg — 
l/\/2/ a - While local strings can be accurately simulated by using the Nambu-Goto action, 



it is difficult to simulate the dynamics of global strings using string-based actions |*j and 



^ While global strings can be modeled by the Kalb-Ramond action, it does not fit for numerical analysis. 
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evolution of the string network has been poorly understood. A crucial feature that makes 
global strings different from local strings is that a long-range force mediated by massless 
NG boson operates between two global strings. 

In general, a scaling solution is characterized by the scaling parameter £, which we 
define as 

£ = E^I^t 2 . (2) 

/^string 

£ is nothing but the average number of infinite strings in a Hubble volume (strictly speak- 
ing, in a box with a volume t 3 ) and should be a constant of order of unity as long as a 
scaling solution is realized. For local strings we find £ ~ 13 from the simulations based on 
the Nambu-Goto action [H] and £ ~ 6 from field-theoretic simulations [15] in the radiation 



dominated universecfl Whereas £ can be much smaller for global strings [T714T9] . 

In addition, the energy stored in global strings is released in a different way from local 
strings. NG bosons dominantly carry away the energy in global strings, while gravitational 
waves play the dominant role for local strings. There has been a controversy about the 
energy spectrum of axions radiated from axionic strings. Davis, Shellard and co-workers 
insist that the spectrum has a sharp peak at the horizon scale [201121] . On the other hand, 
Sikivie and co-workers claim that it is proportional to the inverse momentum [22123] . After 
the QCD phase transition, radiated axions finally become CDM and their energy density 
is proportional to the number density. The energy spectrum of axions is of particular 
importance since the number density of radiated axions is determined by the spectrum. 

A solution for the controversy was given by Yamaguchi, Kawasaki and Yokoyama [T7] 
(hereafter YKY99). In their analysis, a field theoretic simulation of axionic string was 
performed, which is a first-principles calculation and least contaminated by theoretical 
uncertainties. Their analysis showed that the energy spectrum of axions from strings 
are sharply peaked at a momentum around inverse of the horizon scale and suppressed 
exponentially at higher momenta. Therefore the observed shape of the spectrum is in good 
agreement with the insistence by Davis, Shellard and co-workers. In addition, YKY99 
showed that a scaling solution is also realized for global strings and find a much smaller 



scaling parameter £ ~ 1.00 ± 0.08 compared with local stringsP 3 ! This illuminates the 
quantitative difference in dynamics of global strings from local strings. 

Our primary purpose of this paper is to update the analysis of YKY99. While we 
perform a field theoretic simulation of axionic strings of largest scale so far, we also de- 
velop several new techniques to improve the accuracy of analysis. One is a new method 
for identification of strings in a simulation box. For field theoretic simulation, this is a 
non-trivial task since field values are known only at discrete spatial points. By checking 
the consistency with the previous identification methods, reliability of our understanding 
of string dynamics would be much improved. In addition, we introduce a pseudo power 
spectrum estimator (PPSE) [211125] to remove the contamination of string cores in estima- 
tion of the energy spectrum of radiated axions. Since strings are highly-energetic objects, 



# 2 Thc reason of such a discrepancy is discussed and hinted in Ref. |16j . 

# 3 In more refined simulations, a slightly lower value of £ has been obtained [19] . 
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reliable estimate of the energy density of free axions is difficult near string cores. There- 
fore we should remove the regions near strings in estimation of the energy spectrum. In 
YKY99, this was done by dividing the simulation box into eight sub-boxes with the same 
size and estimating the spectrum using only selected sub-boxes with no strings. In order 
to take a larger simulation box and increase the statistics, we need a more effective way 
to remove string cores. This can be done by adopting PPSE, which is often used in power 
spectrum estimation in CMB data analysis to remove observed regions contaminated by 
foregrounds including galactic emission and point sources |26j . 

The organization of the paper is as follows. In Section [2} we first give details of our 
model and the setup adopted in our field theoretic simulation. Then we describe our 
analysis method in Section [31 Particularly, we focus on methods adopted in identification 
of strings and estimation of the energy spectrum. Results of our analysis are presented in 
Section HI where we mainly discuss the scaling property of axionic strings and the energy 
spectrum of radiated axions from strings. Comparison with the previous result of YKY99 
is also discussed here. A constraint on f a and implications for the early Universe are 
discussed in Section [51 The final section is devoted to summary and discussions. 

In this paper, we assume flat Friedmann-Robertson- Walker Universe, with a metric 

ds 2 = -dt 2 + R{t) 2 5 ij dx i dx i , (3) 

where R(t) is the scale factor. We also use a conformal time r = f dt'/R(tf). A dot 
represents derivative respect to the proper time, i.e. ' = d/dt. 



2 Models and setup of field theoretic simulation 

We simulate dynamics of a complex PQ scalar field $(af, t) with a Lagrangian density 

C= |V M $| 2 -K ff [$;T], (4) 



with an effective potential at finite temperature T 



U cff [$;T] = ^(|<&| 



3 1 1 ' 



(5) 



where r\ = f a is the energy scale of PQ symmetry^ 4 I and A is the self-coupling constant. 
The same potential is also adopted in YKY99. When the temperature is high enough 
T > T crit = \/3r), V e s has a minimum at $ = and U(1)pq is restored. At lower 
temperature T < T crit , symmetry breaking minima appear at |$| = rj^l — (T/T crit ) 2 . 
The phase transition is of second order. While we consistently adopt a particular set of 



# 4 In general, f a = r/fN^yj where -/Vdw is the number of degenerate vacua after QCD phase transition. 
Throughout this paper, we assume ./Vdw = 1, so that axionic domain walls quickly disappear after the 
QCD phase transition (otherwise, it is cosmologically disastrous). 
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model parameters t] = 1.22 x 10 16 GeV, A = 1 in our simulation, different choices of these 
parameters do not lead to any qualitative differences in physical consequences. 

We simulate the evolution of <&(x, t) from the initial time t- m i = 0.25t cr i t to the end time 
t cn d = 25t cr it, denoting the time of the PQ phase transition by t cr i t , i-e. T(t crit ) = T crit . 
At the initial time ti n ;, we generate an initial condition for <&(x,t) and t). At high 
temperature T > T crit , $ is approximately in thermal equilibrium, so that $ and $ can be 
regarded as Gaussian random variables. By decomposing $ into the real and imaginary 
parts, i.e. $ = (0i + i(f)2)/v2, correlation functions of <pi and </>2 are given by 

/d 3 k 1 1 

e -^(g-g) (6) 
(27r) d w(jfc,i) e"^*)/^*) - 1 

*)*0 6 (*', t)> = 5 a6 / 7f^e-^^) ^ , (7) 

7 (271") d e u(k,t)/T{t) _ I 

(<j> a (2,t)fo(g,t)) = 0, (8) 

where oj(k,t) = ^ k 2 / R{t) 2 + m{t) 2 is the energy of a Fourier mode with wave number fc, 
with the mass of $ being denoted as m(t), i.e. m{t) 2 = d 2 V cS /d$*d$\<s, =0 = X(T(t) 2 /3 - 
rf). The subscripts a and b can be either 1 or 2. We also note that in Eqs. (jED-flU), we omit 
the contribution of vacuum fluctuations, which would not affect the classical dynamics of 
$. 

The equation of motion for $ is given by 



where H(t) = R(t)/R(t) is the Hubble rate. In our simulation, we integrate Eq. fl9]) using 
the second order leapfrog scheme with a constant conformal time-step At = 2 x 10 _3 r cr i t . 
The background Universe is dominated by radiation. There, the Friedmann equation is 
given by 

= ^Z^ T 4 (10) 

where g* is the number of relativistic degrees of freedom (d.o.f.). In our simulation, we 
take it constant, g* = 1000, following YKY99. 



#5 



Our lattice simulation has the number of grids = 512 3 , which is larger than any 
other previous simulations. We impose a periodic boundary condition and the physical 
size of our simulation box at the end time t en d is taken to 2(r cnd — r crit )/r cnd = 1.6 times 
the horizon scale 1/ H(t cn &). With this choice of the box size, massless axions, which begin 
to be emitted at t cr ; t , travel for a distance less than half the size of our simulation box 
until t cnc j. Therefore, we can avoid axions from overlapping around the box and causing 



# 5 While this value is much larger than that predicted in Standard Model (and most of its extensions), 
it does not affect our main results thanks to the scaling properties of axionic strings. 
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unphysical boundary effects via the long range force mediated by them. In this setup, the 
physical size of lattice spacing at the end time t end turns out to be R(t en( {)Ax = 1.4c? string . 

Finally we comment on the dynamical range of our simulation. On one hand, our simu- 
lation box should be larger than the horizon, whose comoving size scales as l/R(t)H(t) oc 
R(t), in order to avoid boundary effects. On the other hand, the lattice spacing should be 
sufficiently smaller than the string width, whose comoving size scales as 1/R(t); otherwise, 
strings cannot be resolved in the simulation box. As R(t) becomes large, it turns out to 
be impossible to satisfy both of these two demands due to the finiteness of N gr ^- These 
determine the dynamical range. 

3 Analysis method 

Our field-theoretic simulation is a first-principles calculation and free from theoretical un- 
certainties. Several difficulties, however, arise in extracting the physically relevant quan- 
tities from the given data of $ and $ at discrete points. As discussed in Introduction, 
identification of strings and estimation of energy spectrum of axions are of primary im- 
portance. So we describe the details of these matters in this section. 

3.1 Identification of strings 

Several methods for string identification have been developed so far. In Refs. [IT] US], 
a lattice was identified as a part of a string based on the value of the potential energy 
there. Counting number of lattices thus identified as containing strings, they estimated the 
scaling parameter. Although the overall features were traced reasonably well, this method 
had some problems in that string segments were occasionally found disconnected and it 
might overestimate the scaling parameter because the number of the lattices penetrated 
by a string was used rather than the length of the string itself. 

Later, two of the present authors (MY and JY) developed a new method of string iden- 
tification [19], which uses the phase of the fields. In this method, first a quadrate, which a 
string penetrates, is identified by dividing the phase into three zones with unequal spans# 6 
and monitoring the phase rotation just as in the Vachaspati-Vilenkin algorithm [27]. Then, 
the position of the string in each quadrate is determined by the zeros of the two real fields, 
0i and <p2- These complicated identification scheme guarantees the connectedness of strings 
and makes it possible to calculate the string velocity, the intercommutation rate, and the 
NG boson emission rate as well as the scaling parameter more precisely. Thereby it enables 
us to analyze the evolution of the string network from the Lagrangian viewpoint to clarify 
its fundamental characters. 

one divided the phase into three equal zones, one would occasionally identify a quadrate as con- 
taining a string even if either 0i or </> 2 takes the same sign at its four corners. So the authors of [19] had 
to divide the phase into uneven zones. Nevertheless it did not cause any artifacts in the final results. 



5 



For our practical purposes here, however, we do not need to adopt such a complicated 
and time-consuming procedure albeit its accuracy, because the final error is dominated 
by the finiteness of the simulation box rather than the identification scheme we use. We 
therefore adopt a new and much more efficient method in which we do not need to divide 
the phase into fixed uneven zones. 

Let us start by considering a small quadrate whose four vertices are the neighboring 
grids in the simulation box (See the left panel of Figured]). Since the field values t) 
at these vertices (A, B, C and D in the figure) are known from simulation, we can map 
these vertices into the field space (Right panel of the figure). While there are several 
ways to take a region of phase which contains the images of these four vertices, we denote 
the range of minimal one by A6, which is explicitly shown in the right panel of Figure 
CD As shown in the figure, if a string penetrates the quadrate, A9 should be larger than 
7r. This can be easily understood by drawing a line on the quadrate running through the 
penetration point of the string. When the quadrate is small enough for the continuous 
change of the phase to be regarded as isotropic around the string, the phases at opposite 
sides of an arbitrary line intersecting the penetration point (for example, lines of <fii = 
and 02 = shown in the left panel of Figure [I]) differ by ir. Whatever line we draw, 
each region of the quadrate divided by the line contains at least one vertex. Therefore, 
the minimal phase difference cannot be smaller than ir. As long as the continuous phase 
change around strings is isotropic, the opposite is also true, i.e. if AO > ir, the quadrate is 
penetrated by a string. Therefore we can use AO as a criterion for identification of strings. 
One interesting point is that our method is also applicable for any convex polygons other 
than quadrate. In addition, our method is invariant under the global rotation of the phase 
of which is a desirable feature in identification of global strings. 

Of course there are loopholes in our new method. Our method assumes the isotropy of 
the continuous phase change around strings. This is valid when a string is nearly straight 
and there are no other strings nearby. However, this assumption breaks down around re- 
gions where strings are colliding with others and/or strings have small curvature radii. In 
rare cases, two or more strings can penetrate one single quadrate in such regions. There- 
fore, some strings are not connected completely in this method. However, fortunately, 
the fraction of such regions is quite small; it is at most 1% when the system of strings 
relaxes into the scaling regime. Moreover, at around these regions, phase of $ changes 
very violently if not isotropically. Therefore the criterion AO > ir is still useful to identify 
strings, although it to some extent over- or underestimate penetration of strings. In the 
end, we can safely say that the fractions of both omission and commission are less than 
1% in our method. 

When a quadrate is expected to be penetrated by strings, then we determine penetra- 
tion points of strings in it. Our determination method of the penetration points is basically 
the same as in [19]: we first compute the points with <pi = or <p2 = on the boundaries 
of a quadrate using linear interpolation of $; then we identify the intersection of lines that 
connect two points with either <pi = or <f) 2 = as the penetration point (See the left 
panel of Figured]). It sometimes happens that two or more strings are penetrating a single 
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Figure 1: Schematic view of our method for string identification. (Left) Shown is a 
quadrate in real space penetrated by a string (green point). The loci of <pi = and 
4>2 = (red lines) intersect with each other at the string. Note that the string does not 
necessarily penetrate the quadrate perpendicularly. When the phase of $ continuously 
changes around the string isotropically, phases at the opposite sides of an arbitrary line on 
the quadrate intersecting the penetration point differ by it. Greek numbers (I), (II), (III) 
and (IV) in the left panel shows the quadrants of (0i, 2 ) in these regions. (Right) Shown 
is the mapping of four vertices in the field space. The minimal phase range containing the 
images of these four vertices is indicated with a red arrow. By observing whether the phase 
difference A6 in this minimal range is larger than n, we identify quadrates penetrated by 
strings. 



quadrate. In such cases, there are four (not two) points with either <pi = or <p 2 = 
on the boundaries of the quadrate. We compute the sign of <fri or 2 at the center of the 
quadrate by averaging over those at four vertices, and then determine two pairs of points 
each of which we connect so that the sign of 0i or 02 at the center is respected. Sometimes 
intersections of <f>\ — and 2 = are found lying outside the quadrates, but more than 
99 % of penetration points are found lying inside each quadrate. In conjunction with the 
omission and commission rates in identification of penetrated quadrates, we estimate our 
method can determine the positions of strings with at least 99 % accuracy. This is highly 
sufficient to analyze the dynamics of strings and estimate the energy spectrum of radiated 
axions with masking of strings. 
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3.2 Estimation of the energy spectrum of radiated axions 
Energy spectrum of axions 

From the simulated data of $ and $, we obtain the time derivative of the axion field 



a(x, t) = Im 



'ir, 



and its Fourier component is given by 



a(k,t) = I d 3 xe ik *d(x,t) 



121 



where k is a comoving wave number. 

Assuming the statistical isotropy and homogeneity, the two-point correlation function 
in Fourier space can be represented in terms of the power spectrum P(k), 



±(a(k,t)*a(k',t)) = ^5®(k-k')P(k,t), 



(13) 



where the brackets (■) represent an ensemble average, i.e. an average over infinite re- 
alizations. While homogeneity and isotropy is in reality broken by the finiteness of the 
simulation volume, they are approximately valid as long as the wave number k is not 
very small (k < 1/L) nor very large (k > 1/Ax). The mean kinetic energy of the axion 
p(t) is nothing but the ensemble average of energy density p(x,t) = |d(x, t) 2 , and can be 
rewritten in terms of P(k), 



Pit) 



-d(x,t) 2 
d 3 k f d 3 k 



(2tt) 3 



(27T) 



(d(jfe, t)*d(k',t))e^'-^ 



(14) 



This shows that the power spectrum P(/c) defined in Eq. ([TBI) is nothing but the energy 
spectrum of axion field. 



PPSE of the energy spectrum of radiated axions 

What we concern is the energy spectrum Pf ree (/c) of free axions radiated from strings. 
However, d is also induced by moving strings. Separation of contribution from free axions 
and contamination from string motion is quite difficult. Moreover, the energy density 
associated with moving strings is so large that it can easily dominate over that. Therefore 
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it is essential to remove regions around strings in estimation of Pf ree (A;). In the present 
work we adopt PPSE to deal with this issue. 

Let us start from the time derivative of the axion field a(x,t) at a fixed time t. At 
most points in the simulation box, a represents the time derivative of free axions, which 
we denote by df ree . However, near a string, a is also induced by motion of strings, so that 

d(x,t) = dfree + (contamination from strings). (15) 

In the same way as in Eq. ( 1T5j) . Pf ree (k) is defined by 

l -{a lrcc (k,tya, me (k\t)) = ^^5^\k-k')P lrcc {k,t). (16) 

For a while, we drop the argument t because the energy spectrum is calculated at each 
fixed time. 

Fortunately, the contamination is localized around string cores. Therefore we can mask 
the string contamination by adopting a window function 

W(x)={° ^ String8) . (17) 
I 1 (elsewhere) 

If the mask covers regions large enough so that the string contamination is removed to a 
negligible level, masked d has contribution only from free axions. Thus we obtain 

~a(x) = W(x)a(x) = W{x)d {icc {x). (18) 
In the Fourier space, a(k) is the convolution of W(k) and df ree (k), 

~d(k) = J j^W(k - k')d(k') = J j^W(k - k')a ime (k'). (19) 

We can straightforwardly calculate the 'masked' energy spectrum in a given simulation 
box, 

k 2 f dQh 1 ~ - 2 



p{k) s v J 172 4W • (20) 

where V is the comoving volume of the simulation box and the integration is performed 
over the angular direction of k with the solid- angle element dVt^. However, as shown in 
Appendix lA~t the masked spectrum is biased, i.e. (P(k)) ^ P{ rcc (k). This is because the 
wave number k to which d{ rcc (k') contributes can be different from the original one k' ^ k 
due to the mode-mixing induced by the window function W(k — k'). Such a mode mixing 
can be corrected using PPSE. We define a PPSE of Pf ree (/c), 

k 2 f r\k' 

P{k) = V / f^ M_1 (MW'), (21) 
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where M 1 



(k, k') is defined so that it satisfies 



/ 



k' 2 dk 
2vr 2 



M- 1 (fc,fc')M(fc / ,fc") 



2tt 2 



8{k-k"), 



(22) 



(23) 



Then -P(fc) is unbiased, i.e. (P(k)) = Pf Tee (k). The proof is given in Appendix lAl 



Pipeline and validity check 

Figure [2] shows our pipeline. For each single realization of the field theoretic simulation, a 
configuration of a(x) is given at first. Then we identify strings in the simulation box and 



Now we check the validity of our method. For this purpose, we performed a field 
theoretic simulation with lower resolution than that described in Section [2j Here we take 
the lower number of grids (iVgrfd = 256 3 ), decreased dynamical range i en d = 16t CT it and 
smaller box size, 1.0/if (t en d). Other parameters are the same as in Section[2J In estimation 
of the spectrum, we first divide the whole simulation box into eight sub-boxes of the same 
size. Then in each sub-box, we calculated two different spectra estimated with or without 
masking. We note that while the energy spectrum might be affected by boundary effects, 
this does not cause any problem; what we concern here is verification of the method, 
not the physical consequences. In Figure [3], we plot three different spectra estimated at 
t = 16t crit obtained by averaging over 20 independent realizations. The spectrum in red 
crosses corresponds to our PPSE P(k) with masking. The spectra in green boxes and blue 
diamonds correspond to energy spectra calculated without masking of strings. While all 
sub-boxes are used in calculation of the spectrum in green boxes, only sub-boxes found 
without strings in them are selectively used for one with blue diamonds. Therefore the 
blue spectrum is exactly what is calculated in YKY99. The co moving wave number k and 
estimated energy spectra P(k) are normalized in r~ it and r~f t , respectively. 

First of all, we see excellent agreement between our PPSE P{k) (red crosses) and 
the spectrum obtained only from sub-boxes without strings (blue diamonds) in Figure [3j 
Therefore our estimation method of the energy spectrum of free axion can be regarded 
as working quite well. The differences are less than 10 % at /cr crit < 100 and at most 
25 % at /cT cr it < 200. We regard this discrepancy as systematic errors in estimation of 
Pf ree (k). We also observe the energy spectrum estimated without masking (green boxes) 
has significant contribution from string cores at large k. This is because the phase of PQ 
field at some point fixed in a comoving frame rapidly changes when a string passes nearby. 
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Figure 2: Schematic overview of our pipeline. 



However, such a string does not radiate axions and free axions are not responsible for the 
observed large a. Since such large a induced by motion of strings is concentrated around 
strings, it gives power at large k. Since the contribution of string cores is also significant 
at moderate scales /cr cr ; t < 100, it is essential to remove the contamination from string 
cores to estimate the spectrum of energy released by axion emission, which we will discuss 
in Section 14.21 

As a final remark of verification check, we note that at small scales /cr crit > 150, even 
red and blue spectrum differ from the true spectrum which can be calculated from whole 
the simulation box. This is because we performed discrete Fourier transformation in each 
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Figure 3: Validity check of our estimation method using PPSE. Three different spectra 
are plotted (See text for details). Only statistical errors are shown; bars corresponds to 
the square root of the diagonal components of covariance matrices. 



sub-boxes assuming a periodic boundary condition. Since the boundary condition does 
not hold actually within each sub-box, discontinuities in a(x) arise at boundaries and this 
gives additional power to the energy spectrum at large kW\ True energy spectra should 
have to some extent lower amplitudes at large k (See Figure EJ). 



4 Results 

4.1 Scaling property of axionic strings 

Figure H] is a visualization of one realization in our simulation. There, red points are the 
positions of axionic strings determined by the method explained in the previous section. 
The points clearly form line-shaped objects regarded as axionic strings, which shows good 
resolution in our simulation as well as excellent string identification. 

# 7 This is the very reason that the energy spectrum shown in Figure 2 of YKY99 shows deviation from 
exponential behavior at large wave numbers, which is not mentioned in YKY99. 
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Figure 4: Visualization of one realization from our field theoretic simulation. Red line 
corresponds to axionic string identified by our method discussed in Section 13.11 r in the 
top right of each panel is the conformal time of each time slice, which can be translated 
into the proper time t via a relation in radiation domination t/t cr - lt = r 2 /r c 2 rit . The spatial 
scale shows a co moving length in unit of the horizon size at t C nd = 25t cr i t . 



By linearly connecting the penetration points of strings in neighboring quadrates, we 
can estimate the total length of strings in the simulation box. Recalling that p s tring//4tring 
in Eq.(j2]) corresponds to the mean physical length of strings in unit physical volume, we 
can compute the scaling parameter £ from the total length of strings. In Figure [5j we 
plotted the time evolution of £ which is obtained by averaging over 20 realizations with 
the setup explained in Section [2J We first see that £ stays constant for t > 10t cr i t . This 
shows that the system of axionic strings relaxes into the scaling regime. At t — 25t cr ^, we 
obtain 

£ = 0.87 ±0.14. (24) 
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Figure 5: Time evolution of the scaling parameter £ obtained by averaging over 20 real- 
izations. Note that data points are not homogeneously placed in t, but in r oc \fi. 



Our result shows good agreement with [T5], which gives £ ~ 0.8. While YKY99 gives a 
slightly higher value £ = 1.00 ± 0.08, the difference is not significant at all. 

It might be curious that the error in £ grows as t increases in Figure [5j This is due to 
a statistical reason; since the horizon scale becomes larger as time advances, the number 
of independent horizon volumes in the simulation box effectively becomes smaller at later 
times. Therefore, variation in £ among realizations becomes larger as t increases, which 
increases the errors in £ averaged over realizations. Moreover, our estimate of £ in Eq. 
([24"]) can be regarded as rather conservative one. 



4.2 Net energy spectrum of radiated axions 

Figure [6] shows the energy spectra of radiated axions at t± = 12.25t crit (left) and t 2 = 25t crit 
(right), that are estimated from the 20 realizations used in Section [4.11 The amplitude 
of energy spectrum at t 2 is about (ti/t 2 ) 2 — 0.24 times that at t±. This is because the 
energy density of free axions scales as R~ 4 (t), without emission or absorption. We see a 
clear exponential behavior at large k after the removal of the contamination from strings. 
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Figure 6: Energy spectra of radiated axions at t\ = 12.25t crit (left) and t 2 = 25t crit (right). 
Green bars correspond to statistical errors alone; red error bars also count systematic 
errors suggested in Section I3T21 as well as statistical ones. 



As we see in Section H~Tj the system of axionic strings are already in the scaling regime. 
Most of axions at this epoch are however emitted before the settlement into the regime. In 
order to extract the energy spectrum of axions radiated during the scaling regime, we need 
to differentiate the energy spectra at different times. We define the differential spectrum 
of radiated axion between t\ and t 2 , 

AP ilcc (k;t u t 2 ) EE R\t 2 )P ilcc (k,t 2 ) - i^tOPfceeCMi). (25) 

If there are no emission nor absorption of axions, the energy density of axion scales as 
R~ 4 (t). Therefore, APf ree (fc, t) is the net energy spectrum of axions radiated from strings. 

In Figure we plotted the differential spectrum between t\ and t 2 . We observe that 
the net energy spectrum is sharply peaked at around the horizon, which corresponds to 
k = 3.6r CTi ^ at t — 1\ (k = 2.5r~^ at t = t 2 ), and its amplitude is exponentially suppressed 
toward higher k. This is consistent with YKY99. 

To calculate the number density of axions created from strings, the important quantity 
is the mean reciprocal comoving momentum of radiated axion defined by 



fc-i(i) = v.; (26) 



By fitting the ratio of R{i)k 1 (t) to the horizon scale £/27r, we obtain 



R(t)k- l (t) i . . 
—— ee 6 = 0.23 i 0.02. 27 

This is consistent with e _1 = 0.25 ± 0.18 obtained by YKY99. However, the uncertainty 
is reduced by an order of magnitude owing to the increase in statistics, which is brought 
about by use of PPSE. 
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Figure 7: Differential energy spectrum of radiated axions between t\ = 12.25t cr it and 
£2 = 25t cr i t . Errors shown are estimated from the quadrature of those shown in Figure El 
Green (red) bars correspond to statistical errors alone (statistical and systematic errors). 
Note that the scale in y-axis is arbitrary since the scale factor R(t cri t)/Ro is not well- 
defined. 



5 Constraint on the axion decay constant 

Using the result of the previous section, we can derive a constraint on f a . After the PQ 
phase transition, the system of axionic strings reaches the scaling regime and continuously 
emit massless axions. After the QCD phase transition, the axionic string-wall system is 
generated. After t w , when the tension of axionic domain walls a" wa u starts to dominate the 
system, i.e. cr wall (i w ) = fJ, s trmg(tw)/tw, the wall-string system quickly disappears. Emission 
of axions from strings continues until t w , and the largest fraction of the axions are produced 
just before t w . After t w , axions becomes non-relativistic due to the finite mass. Today, 
axions exist as CDM in the Universe. 

In Appendix |Bl we give a detailed derivation of the density parameter of CDM axions 
radiated from axionic strings, adopting the axion mass at finite temperature from recent 
studies j7J[33]. By substituting Eqs. ( 124)) and (127)1 obtained from our simulation into Eq. 
( 150)) . we obtain 

ftaxion/* 2 = (1.66 ± 0.25)7 ( ^) ~°' 31 ( 1 ( — ^ 1 9 , (28) 

where 7, g* w and A are the dilution factor, the number of relativistic d.o.f. at t w , and the 
scale of the QCD phase transition, respectively (See Appendix IB")) . 
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^axion^ 2 should be smaller than the observed f2cDM^ 2 - Recent cosmological observa- 
tions [28] give ficDM^ 2 = 0.11 (assuming flat power-law ACDM model). Therefore we can 
translate Eq. ( J28l) into constraint on f a . Assuming no entropy dilution (7 = 1), we obtain 
f a < 1.3 x 10 11 GeV at 2 a level. By taking account of uncertainties in the QCD phase 
transition, a conservative constraint would be 

f a < 3 x 10 11 GeV. (29) 

Our constraint is a few times tighter than one given in YKY99. This is partially because 
of the precise Ocdm^ 2 obtained from recent observations, which was not available at the 
time of YKY99, and also of the axion mass at finite temperature from [H[33], which is 
slightly smaller (hence larger t w ) than that adopted in YKY99. 

We also note that while the error in £ of Eq. ( 1241) can be smaller if it is estimated at 
earlier times t < t en d, this does not affect our final consequences much; the error in the 
energy density of axion CDM in Eq. (|28|) can decrease at most 40% and the upper bound 
on f a does not change. 

6 Summary and discussions 

Using field theoretic simulation, we have investigated the dynamics of strings and the 
energy spectrum of axion radiation from strings, which was previously done by YKY99. 
Our lattice simulation has a higher resolution than that of any other previous studies. 
We have also introduced several new techniques in order to improve accuracies of the 
previous analysis. We have developed a new identification method of strings in simulation 
box, which uses the minimum phase difference among neighboring grid points in quadrates. 
The estimated scaling parameter £ = 0.87±0.14 shows good agreement with other previous 
studies. The consistency among completely different identification methods would suggest 
that these results are conclusively robust. 

In estimation of the energy spectrum of axions radiated from strings, we have also 
introduced a new method to obtain larger statistics. We have verified our new method 
using PPSE, by checking consistency with the spectrum calculated from sub-boxes in the 
simulation box found without strings. The differential energy spectrum obtained from 
our lattice simulation peaks sharply around a horizon scales and damps exponentially 
toward higher wave numbers. This shows good agreement with YKY99, and support the 
discussion of Davis and Shellard [20]. The ratio of the horizon scale and the mean energy 
momentum of radiated axion is e _1 = 0.23 ± 0.02. Using these results, we have obtained 
a constraint on the axion decay constant / a <3x 10 11 GeV. 

Above constraint on f a counts only axions radiated from strings. However, axions 
are also radiated from axionic domain walls. This gives a constraint f a < 2.5 x 10 11 
GeV [2"3" jl2"9"H3"2"] . which is also severer than in the original papers due to the more precise 
observationally inferred value of ^cdm^ 2 - Moreover, the oscillation of zero mode of axion 
field also contributes to the energy density of CDM axions in the present Universe, which 
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gives f a < (2.8 ±2) x 10 11 GeV [TJ. Combination of these constraints would give f a < 10 
GeV, as long as there occurs no entropy dilutions after the QCD phase transition. 
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A Unbiasedness of pseudo power spectrum estimator 

Here we prove that the pseudo power spectrum estimator P(k) given in Eq. (12 ip is an 
unbiased estimator of the energy spectrum P lrcc {k) defined in Eq. (TIB"]) . Since a is a 
convolution of W and af ree , the masked spectrum in Eq. (|20|) can be rewritten as 

m = j^J^j*£L\w(i- k'yw(k - ?wfrwf). (30) 

Using Eq. (116j) . the ensemble average of the masked spectrum can be written in terms of 

-Pfree(^); 

/6„» [<M k e [ d 3 k' 1 r r,. 2 

{m} = J 17V J Wfk* w(k ~ k } P '"" (k ' 

/rlk' 
^M(k,k')P ilcc (k'), (31) 

where we have used Eq. (12"3"|) in the second equality. Unless there are no strings and 
W{k) = (2ir) 3 5^(k), (P(k)) ^ Pf ree (fc). Therefore P(k) is in general a biased estimator 
of Pf Iee (k). 

On the other hand, the ensemble average of PPSE can also be calculated. From Eq. (|21J) 
we obtain 

(P(k)) = y J I^M-\k,k')(P(k>)) 

= yj ^M-\k,k')Vk' 2 J ^ 2 M(k',k")P iiee (k") 

= Piree{k), (32) 

where in the second and third equality, we used Eqs, (|3~T|) and ff22l) . respectively. Therefore 
P(k) is an unbiased estimator of Pf ree (fc). 
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B Energy density of axion CDM from strings 

The energy density of axions radiated from strings are to be calculated. When a network 
of axionic string is in a scaling regime, mean energy density of strings is given by 



PstringO) = ^2vr/ a 2 ln 



tl>/l 



d 



string 



(33) 



where in the right hand side, a factor £/t 2 represents the mean density of physical length 
of strings, and the rest 2nf^ ln(i/\/£d s tring) is the line energy density of strings. Strings 
release their energy mainly by emitting massless axions. Therefore the evolution equations 
for the system of strings and radiated axions are given by 



^Pstring (t) 
Jt 

^Paxion(^) 

dt 



-2H (t)p s tring(£) — 

-4#(t)p axion (t) + 



^Pstring (t) 
Jt 

G^Pstring (t) 

dt 



emission 



(34) 
(35) 



J emission 



Using Eq. (]33p . the energy loss rate of strings via axion emission is obtained from Eq. 



^Pstring(^) 
Jt 



emission 



hi 



string 



We define 



N(t) = R(t) 3 n axion (t), 
E(t) = R(t) 4 p axion (t). 



(36) 



(37) 
(38) 



N(t) is the number of axions in a unit comoving volume. Both N(t) and E(t) are constant 
for massless axions if there are no production or absorption of axions. Then, Eq. (|35|) can 
be rewritten in terms of E(t). We obtain 



dE(t) 
dt 



R(tY 



^Pstring(^) 
Jt 



(39) 



J emission 



The number of axions radiated from strings can be given by 

N(t) 

where we have used 



v ; df 



(40) 



dN(t)/dt _ Tt 
dE{t)/dt ' 



dt 



[R(t) 4 f |^Pf r ec(M 
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k-i(t). 



(41) 



t* is the time when scaling solution is realized. Assuming is proportional to the 

horizon scale 

1 t 



R{t)k-i(t) = --, 
e 271 



(42) 



we obtain 



N(t) 



dt' —R(t') A 2nfX— 



In 



d, 



string 



hi 



d, 



string 



(43) 



where in the last line, we have used that t oc R(t) 2 is a good approximation in the radiation 
domination. Strictly speaking, we need to take into account the change of the degree of 
the freedom of relativistic particles. However, axion emission from strings continues until 
the time of wall formation t w ^> t* and actually is dominated by the contribution at the 
last time t w . Thus, we can safely omit such change of the degree of the freedom. After t w , 
the number of axions in a comoving volume is conserved, so that 

\3 



N(t > t w ) = N(t w ) ~ 2fl 



2 £ R{t % 



1 1, 



In 



(44) 



The number density of axions from strings are ri ax i on (to) = N(t w )/R(to) 3 . 

Regarding the axion mass at finite temperature, we quote the recent result from 
where the Interacting Instanton Liquid Model is adopted. If a simple power-law fit is 
applied, the axion mass is given by 



A 4 (T 
m a (T) = a— - 



Pa VA 



(45) 



where n = 6.68, a = 1.68 x KT 7 and A = 400 MeV. At small T, Eq. (ggj) can be arbitrary 
large, and when it becomes larger than the axion mass at zero temperature 



A 4 

m a {T = 0) 2 = 1.46 x 10~ 3 -^ 

J a 



2 • 



(46) 



we simply set m a (T) = m a (T = 0). This approximation gives good agreement with other 
studies [34"|l3"5] . 

At the time t w , the tension of axionic domain wall <J wa \\{t) begins to dominate the 
axionic string- wall system, i.e. a wa \\(t w ) = fi s trmg{t w )/t w . The wall tension is given by 
Cwaii(^) — 37rm a (t)/ ( 2 . We assume a constant number of relativistic d.o.f. g* around t w , 
i.e. t w = \/2H{t w ). Then from Eq. ( 1431) . t w and the temperature T w at t w are given by 

-2 / r \ 4/(n+4) 



til 



2.2 x 10 



-6 / 9*w 

sec 

70 



-n/2(n+4) 



(■ 



A \ 



T w = 0.64 GeV 



g*w 
To" 



-l/(n+4) 



V400MeV / 
A 



fa 



400MeV/ \10 12 GeV 



fa 



10 12 GeV / 

2/(n+4) 



(47) 
(48) 
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where g* w is g* at t w . When hadronic states up to a mass of 3 GeV are counted [7j, we 
find g* w ~ 70 for T w ~ 1 GeV. 
Entropy conservation yields 

R( T w )^7.2 xW -«,(^)^y 3 , (49) 

where 7 is the dilution factor due to the entropy production after t w . Finally, combining 
Eqs. dS]), (l4"7|) and (I4"9l . we obtain the density parameter of axion radiated from axionic 
strings 

^ = 8 - ? 7 7 ( 70 ) li^OMeV J vT^GeV J ' (50) 

where h is the reduced Hubble constant, i.e. H = 100/i km/sec/Mpc. 
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